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Ab initio techniques based on density functional theory in the projector-augmented-wave imple- 

p^ . mentation are used to calculate the free energy and a range of other thermodynamic properties of 

liquid iron at high pressures and temperatures relevant to the Earth's core. The ab initio free en- 
ergy is obtained by using thermodynamic integration to calculate the change of free energy on going 
from a simple reference system to the ab initio system, with thermal averages computed by ab initio 
(^ ' molecular dynamics simulation. The reference system consists of the inverse-power pair-potential 

t^ , model used in previous work. The liquid-state free energy is combined with the free energy of hexag- 

onal close packed Fe calculated earlier using identical ab initio techniques to obtain the melting curve 
and volume and entropy of melting. Comparisons of the calculated melting properties with experi- 
mental measurement and with other recent ab initio predictions are presented. Experiment-theory 
comparisons are also presented for the pressures at which the solid and liquid Hugoniot curves cross 
jrt I the melting line, and the sound speed and Griineisen parameter along the Hugoniot. Additional 

comparisons are made with a commonly used equation of state for high-pressure/high-temperature 
Fe based on experimental data. 
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The last few years have seen important progress in calculating the thermodynamic properties of condensed matter 

using ab initio techniques based on density-functional theory (DFT) [|l|-Q|. There has been particular attention to 

the thermodynamics of crystals, whose harmonic free energy can be obtained from phonon frequencies computed 

by standard DFT methods |q-0|. The ab initio treatment of liquid-state thermodynamics is also important, and 

^Kf-\ ' thermodynamic integration has been shown to be an effective way of calculating the DFT free energy of liquids W,0,0 . 

f — , These developments have made it possible to treat phase equilibria, including melting properties, by completely ab 

f^ initio methods. We report here DFT free-energy calculations on high-pressure/high-temperature liquid iron, which 

^~~l , we combine with earlier results on the solid ||ll| to obtain the complete melting curve and the variation of the volume 

^^ ■ and entropy of melting along the curve. We also present results for some key thermodynamic properties of the liquid, 

which we compare with data from shock experiments and other sources. The general methods developed here may 

be useful for other problems involving phase equilibria under extreme conditions. A brief report of this work was 

presented earlier Q. 

rj^ . The properties of high-pressure/high-temperature Fe are of great scientific importance because the Earth's core 

r"| ' consists mainly of Fe, with a minor fraction of light impurities |l3-n5[. The melting curve is particularly important, 

O , since it provides one of the very few ways of estimating the temperature at the boundary between the liquid outer core 

O ■ and the solid inner core |l6[. Because of this, strenuous efforts have been made to measure the melting curve [0-E3|, 

^ [ but the extreme pressures and temperatures required {p ~ 330 GPa, T ~ 6000 K) make the experiments very 

•""j ' demanding. Ab initio calculations therefore have a major role to play, and several independent attempts to obtain 

rN I the melting curve using different ab initio strategies have been reported recently ||l^,^,^ . The rather unsatisfactory 

S . agreement between the predictions makes a full presentation of the technical methods all the more important. 

- - ' The calculation of melting properties using ab initio free energies was pioneered by Sugino and Car |l| in their work 

on the melting of Si at ambient pressure. Related methods were subsequently used by de Wijs et al. [Q to study the 

melting of Al. In both cases, thermodynamic integration (see e.g. Ref. [p7|) was used to obtain the ab initio free 

energy from the free energy of a simple reference system, and we follow the same strategy here. The other recent 

calculations J25U2q] on the high-pressure melting of Fe employed ab initio methods in a different way. Free energies 

were not calculated, but instead an empirical parameterised form of the total-energy function was fitted to DFT total 

energies calculated for representative configurations of the solid and liquid. The empirical energy function was then 

used in molecular dynamics simulations of very large systems containing coexisting solid and liquid. 



The detailed DFT techniques used in this work are identical to those used in our work on h.c.p. Fe |T^. In 
particular, we use the generalised gradient approximation (GGA) for exchange-correlation energy, in the form known 
as Perdew-Wang 1991 |22,|2^, which reproduces very accurately a wide range of experimental properties of solid 
iron, as noted in more detail elsewhere |30|-|33|]. We also use the projector-augmented- wave (PAW) implementation 
of DFT |p3|-p5[ , which is an all-electron technique similar to other standard implementations such as full-potential 
augmented plane waves (FLAP W) [p6[ , as well as being closely related to the ultrasoft pseudopotential method [^ . 
We have used the VASP code | |3q , |39| , which is exceptionally stable an efficient for metals, with the implementation 
of an extrapolation of the charge density which increases the efficiency of molecular dynamics simulations by almost 
a factor of two Q . 

The calculation of melting properties demands very high precision for the free energies of the two phases, as 
emphasised elsewhere P,|lli . The required precision is set by the value of the entropy of melting, and one finds that 
in order to calculate the melting temperature to within 100 K the non-cancelling error in the free energies must be 
reduced to ~ 10 meV/atom. The use of identical electronic-structure methods in the two phases is clearly necessary; 
but it is certainly not sufficient, since the detailed free-energy techniques differ in the two phases. In the solid, we 
relied heavily on harmonic calculations, whereas the liquid-state calculations rely on relating the free energy to that 
of a reference liquid. It is therefore essential to reduce the statistical-mechanical errors below the tolerance, and we 
aim to demonstrate that this has been achieved. 

In the next Section, we summarise the technical methods, and Sec. 3 then reports our results for the DFT free 
energy of liquid Fe over a wide range of thermodynamic states. Sec. 4 presents our calculated melting properties, 
which we compare with experimental results and the predictions of other ab initio calculations. Our free-energy 
results have been used to compute a variety of other thermodynamic quantities for the liquid, and we compare these 
in Sees. M and VI with direct shock measurements as well as published extrapolations of other experimental data. In 



the final Section, we give further discussion and a summary of our conclusions. 

II. TECHNIQUES 

The key thermodynamic quantity calculated in this work is the ab initio Helmholtz free energy F, which, with the 
statistical mechanics of the nuclei treated in the classical limit, is: 

i^ = -fcBrin|^^^j^ /"dRi...dRjvexp[-/3f/Ai(Ri,...R7v;T;i)]| , (1) 

where R^ (i = 1, . . . N) are the positions of the N nuclei, A — h/{2TTMkBTy^^ is the thermal wavelength, with M 
the nuclear mass and (3 = l/k^T. The quantity [/ai(Ri, . . .RAr;Tei) is the DFT electronic free energy calculated 
with the A'^ nuclei fixed at positions Ri, . . . Rjv- This is given, following the Mermin formulation of finite-temperature 
DFT |^l[, by C/ai = E — TS, where the DFT energy E is the usual sum of kinetic, electron- nucleus, Hartree 
and exchange-correlation terms, and S is the electronic entropy, given by the independent-electron formula: S — 
— fceTci J2i [fi In /i -f (1 — /i) ln(l — /;)] , with fi the thermal (Fermi-Dirac) occupation number of orbital i. In exact 
DFT, the exchange correlation (free) energy E^c has an explicit dependence on Td but we assume here that E^xc has 
its zero-temperature form. It was shown in Ref. |11[| that quantum corrections to the classical approximation to the 
free energy are negligible in the high-temperature solid, and the same will be true a fortiori in the liquid. 

Our earlier work |Q should be consulted for technical details of the PAW implementation of DFT that we use. We 
note here that under Earth's core conditions it is not accurate enough to neglect the response of the 3p electrons, and 
in principle these should be explicitly included in the valence set along with the 3d and the 4s electrons. However, as 
shown earlier l32,p3| , the computational effort of including 3p electrons explicitly can be avoided with almost no loss 
of precision by mimicking the effect of the 3p electrons by an effective pair potential. The procedure for constructing 
this pair potential was described in Ref. ||33| . The pair potential used here is exactly the same as we used in our 
thermodynamic calculations on the h.c.p. solid, so that a good cancellation of any residual errors is expected. To 
avoid any possible doubt on this score, we have also done spot checks on the effect of including both 3s and 3p 



electrons explicitly in the valence set, as reported in Sec. IV. The outermost core radius in our PAW calculations is 
1.16 A. At Earth's core pressures, the atoms in both liquid and solid come closer than the diameter of the ionic cores, 
which therefore overlap. We will show in Sec. 1^ that this too has only a very small effect on the free energies. 

The present calculations, like those reported earlier on the h.c.p. solid, make use of 'thermodynamic integration' p7| , 
which is a completely general technique for determining the difference of free energies Fi — Fq of two systems whose 
total-energy functions are Ui and Uq- The basic idea is that Fi~Fq represents the reversible work done on continuously 
and isothermally switching the energy function from Uq to Ui. To do this switching, a continuously variable energy 
function U\ is defined as: 



Ux = il- \)Uo + XUi , (2) 

so that the energy goes from Uq to C/i as A goes from to 1. In classical statistical mechanics, the work done in an 
infinitesimal change dX is: 

dF = {dUx/dX)xd\ = {Ui - Uo)xdX , (3) 

where {■)x represents the thermal average evaluated for the system governed by U\. It follows that: 

Fi^Fo= f dX (C/i - C/o)a . (4) 

Jo 

We use this technique to calculate the ab initio free energy Fai of liquid Fe by identifying Ui as the ab initio total 
energy function [/ai(Ri, . . .Rat) and Uq as the total energy C/rcf(R-i, • ■ ■'R-n) of a simple model reference system, 
whose free energy Fref can be calculated. Then Fai is given by: 

FaI - ^rcf + f dX {UaI - C/rcf)A ■ (5) 



In practice, we calculate {Uai ~ t/rof)A for a suitable set of A values, and perform the integration numerically. The 
average {Uai — C^ref)A is evaluated at each A using constant-temperature ab initio molecular dynamics, with the time 
evolution generated by the total-energy function U\ = {1 — X)Urcf + XUai- 

As explained in more detail elsewhere |0] , the computational effort needed to perform the thermodynamic integra- 
tion is greatly reduced if the fluctuations of the energy difference AC/ = Uai — Urd are small, for two reasons. First, 
the amount of sampling needed to calculated {AU)\ to a given precision is reduced; second, the variation of (AC/) a 
as A goes from to 1 is reduced. In fact, if the fluctuations 5AU = AC/ — (AC/)ai are small enough, one can neglect 
this variation and approximate Fai — F^f -I- (AC/)ai, with the average taken in the ab-initio ensemble. If this is not 
good enough, the next approximation is readily shown to be: 

Fai ~ Fref + (AC/)ai + -^{{SAU)^)ai • (6) 

2/eBi 

Our task is therefore to search for a model C/rcf for which the fluctuations Uai ~ t'ref are as small as possible. 

The problem of mimicking the fluctuations of ab initio energy Uai in high-p/high-T liquid Fe using a reference 
system was studied in detail in a recent paper [ p3[ . We showed there that a Urd consisting of a sum of pair potentials: 

C/ref = Uth + Upair , (7) 

in which 
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C^pah--2 2^'^(|R.-R, I), (8) 

can be arranged to mimic the fluctuations of C/ai very precisely, if we choose (l)(r) to be a repulsive inverse-power 
potential (f>{r) — B/r", with suitable values of B and a. In this expression for C/rof , we have included a term C/th which 
depends on thermodynamic state, but does not depend on the positions R^. We define this C/th so as to minimize the 
mean square value of AC/, which requires the following condition: 

(AC/)ai = (C/ai - C/pair - Uth)Ai = . (9) 

It might at first seem surprising that such a simple Ui-d can reproduce Uai accurately, since it does not explicitly 
represent the metallic bonding due to partial filling of the Sd-band, which is absorbed into the term Uth- However, 
the only requirement on U^d is that the fluctuations 6AU should be small for the given liquid state, and we shall 



demonstrate below that this is indeed achieved. We comment further on this question in Sec. VI] . 
The free energy of the reference system can be expressed as: 

Frd = C/th + Fpair , (10) 

where Fpair is the free energy associated with C/pair. An important advantage of our chosen form of U^ef is that Fpair 
depends non-trivially only on a single thermodynamic variable, rather than depending separately on temperature 
T and atomic number density n. Let F^^^^ = Fpair ~ ^pg be the excess free energy of the reference system, i.e. 



the difference between Fpair and the free energy Fpg of the perfect gas at the given T and n. Then the quantity 
/pair — ^pair/-^^BT depends Only on the dimensionless thermodynamic parameter C defined as: 

C = Bn-^'^/kjiT . (11) 

This simphfies the representation of i^rcf , since we can write; 

Fref = f/th + Fpg + 7VfcBr/p-,i,(C) . (12) 

We also expect the representation of f/th to be simple. Note first that since t^th is defined so that {AU)ai — 0, then 
from Eqn (0) we must have: 

C/th = {UaI - f/pair)AI ■ (13) 

Now if the fluctuations of Uai — C^pair are indeed small, then the value of {Uai — f/pair)Ai should be very close to the 
value of Uai ~ C^pair evaluated at zero temperature with all atoms on the sites of the perfect h.c.p. lattice having the 
same density as the liquid. Denoting by f/^j and C/p^ir ^^^ values of the zero-temperature ab initio and pair-potential 
energies for the perfect lattice, and defining U^^^ = [/^i ~ t^pair' '^^ can then write: 

f/th = t/t°h + SUtb , (14) 

where SUth will be a small quantity depending weakly on volume and temperature. The accurate computation and 
representation of C/^i were discussed in Ref. Q , and the accurate computation of U^^^^ is clearly trivial, so that the 
treatment of C/^j^ is straightfo rward . The small difference SUth = {Uai — C/pair)Ai — U^-^^ is evaluated from the AIMD 



simulations described in Sec. [II C 



We conclude this Section by summarising our route to the calculation of the ab initio free energy Fai of the liquid. 
We employ Eqn (g), ignoring the higher-order fluctuation terms. Recalling that (AC/)ai — 0, and combining Eqns ([l^ ) 
and dlj), we have: 

FAi^F,,f + {{5AUf)Ai/2kBT 

= [/°h + SUth + Fpg + 7VfcBT/p-,i,(C) + {{SAU)^)Ai/2kBT . (15) 

We now turn to the calculation of the reduced free energy /pair(C) of th^ reference system and the small quantities 
SUth and {{SAU)'^)ai- We shall also give evidence that with our chosen reference model the higher-order fluctuation 
terms omitted from Eqn dlq) are indeed negligible. 

III. FREE ENERGY OF THE LIQUID 

A. Inverse-power reference system 

The PAW calculations used to validate the inverse-power reference system are those reported in Ref. ^^. They 
consist of a set of AIMD simulations performed at 16 thermodynamic states covering the temperature range 3000 — 
8000 K and the pressure range 60 — 390 GPa. All the simulations were performed on a 67-atom system using F-point 
sampling, with a time step of 1 fs. We stress that such a small system with such limited sampling cannot be expected 
to yield very precise results for thermodynamic quantities, and our only purpose here is to demonstrate the adequacy 
of the reference system. At each thermodynamic state, the system was equilibrated using the reference system itself, 
and AIMD data were then accumulated for a time span of 5 ps. 

We showed in Ref. [Q that the inverse-power model, with parameters a = 5.86 and B chosen so that for r = 2 A the 
potential 4i{r) is 1.95 eV, reproduces very closely the ab initio liquid for the state T = 4300 K, p = 10700 kg m^"^. We 
have studied the strength of the SAU fluctuations for all 16 thermodynamic states, using exactly the same reference 
model for all states, and we report in Table S the normalized strength of these fluctuations, which we characterize 

by the quantity a = [[{SAU)'^)ai/N'\ . Two points should be noted: First, a is small, since its typical value of 
100 meV is markedly smaller than the typical thermal energies /cbT (258 meV at the lowest temperature of 3000 K). 
Once a is as small as this, little is gained by further improvement of the reference system. Second, a does not vary 
strongly with thermodynamic state, so that the reference system specified by the values of a and B given above can 
be used for all the thermodynamic states of interest here. 



B. Free energy of reference system 

In order to cover the range of thermodynamic states of hquid Fe that interests us, we need accurate values of the 
excess free energy of the reference system for C values going from 2.5 to 5.0. There have been many studies of the 
thermodynamic properties of inverse-power systems, including one on the free energy of the liquid 1/r^ system [ [42| , 
but since these do not provide what we need we have made our own calculations of /^^f (C) for the l/r^-^^ case. Our 
strategy is to start with standard literature values for the excess free energy of the Lennard- Jones (L J) liquid (we use 
the results reported in Ref. ||43|1 ), and to use thermodynamic integration to go from the LJ system to the inverse-power 
system, so that Uq and Ui in Eqns (^-^) represent the LJ and inverse-power total energies respectively. In doing 
this, our target was to keep technical errors small enough so that the final free energy Fj-cf is correct to better than 
5 meV/atom. 

We note the following technical points. The calculations were done at a standard volume per atom, usually taken to 
be 8.67 A^, with the temperature chosen to give the required value of C- Ewald techniques were used to avoid cutting 
off the inverse-power potential at any distance - we regarded this as essential, since a cut-off would compromise the 
scaling properties of the reference system. The classical molecular dynamics simulations used to compute {Ui — Uo)x 
were done using the constant-temperature technique, with each atom taken to have a mass of 55.86 a.u., and the 
time-step set equal to 1 fs. For each thermodynamic state, we are free to choose any convenient values for the LJ 
parameters e and a. Our criterion for choosing these is that the fluctuations of Ui — Uq should be kept reasonably 
small, but with the proviso that the initial LJ system must be in the liquid state. In many cases, we have checked 
for consistency by using different e and a values. Since we require f^^[{C) in the thermodynamic limit of infinite 
system size, we have made careful checks on size effects. Tests on systems containing up to 499 atoms show that size 
errors in /^^f (C) are less than 1 meV/atom, and this is small enough to ensure that i^rof has a precision of better than 
5 meV. Most of this error arises from the error in the literature values of the LJ free energy. As a further check on 
our techniques, we have done calculations on the 1/r^ system at selected thermodynamic states, and compared with 
the free energy results of Laird and Haymct |4^ . 

We have checked our procedures by repeating most of the calculations using the perfect gas as reference system, so 
as to be free from possible errors in the free energy of the LJ system. For these calculations we used a different form 
for U\, namely 

t/A = (l-A2);7o + A2[/i, (16) 

where Ui is the total energy of the inverse power system and f/o = is that of the perfect gas. Using this functional 
form Eq. ^ becomes 

Fi^Fo= f dX2X{Ui-Uo)x. (17) 



The advantage of using this different functional form for U\ is that the value of the integrand does not need to be 
computed for A = 0, where the dynamics of the system is that of the perfect gas. In this case, since there are no 
forces in the system, there is nothing to prevent the atoms from overlapping, and the potential energy Ui diverges. 
Not computing the integrand at A = only partially solves this problem, since for small values of A the forces on the 
atoms are small, the atoms can come close together, and the potential energy Ui fluctuates violently. However, we 
found that by performing long enough simulations, typically 1 ns, we could calculate the integral with an accuracy 
of ca. 1 meV/atom. These calculations with the perfect-gas reference system give excess free energies of the inverse 
power system that are systematically 5 meV/atom lower than those obtained using the LJ reference system. Our 
belief is that the discrepancy arises from a small systematic error in the free energies given in Ref. Eq| . 

After all these tests, calculations of /j^f (C) were done at a regularly spaced set of ( values at intervals of 0.25, and 
we found that the results could be fitted to the required precision by the following 3rd-degree polynomial: 

3 

/rrf(C)=E^^^^- (18) 

The values of the coefhcients are: cq = 1.981, a = 5.097, C2 = 0.1626, cg = 0.009733. 

C. From reference to full ab initio 

To achieve our target precision of 10 meV/atom in the ab initio free energy Fai of the liquid, two sources of error 
must be studied: system size effects and electronic fc-point sampling. An important point to note is that these errors 



only affect the small terms SUth and {(SAU)'^) ai /'^ksT in Eqn ([l5|), since /,^f (C) refers already to the infinite system, 
and fc-point errors in U^^^ are negligible. We also study the validity of neglecting the higher-order fluctuation terms 
inEqn(|l|). 

We focus first on the quantity 6Uth in Eqn (|l5|). To study size errors in this quantity, we calculated the thermal 
average (C/ai — C^pair}Ai for a range of system sizes. These test calculations were done on systems of up to 241 atoms at 

V = 8.67 A^/atom and T — 4300 K using F-point sampling. The preparation and equilibration of these systems were 
done using the inverse-power reference system. Since the latter so closely mimics the ab initio system for the 67-atom 
cell, it should provide a well equilibrated starting point for ab initio simulation of larger systems. The duration of all 
the ab initio simulations after equilibration was 1 ps. The results of these tests are summarised in Table El, where 
we report the value of 5Uth per atom, i.e. the quantity SUth/N = [{Uai — t/paii)Ai — U^^]/N (see Sec. ^. Since 
U^^/N is independent of the system size, the variation of the reported quantity arises solely from size dependence of 
{Uai — Upa,ii)Ai/N. We see that with ~ 125 atoms SUth/N is converged to better than 5 nieV/atom, and that already 
with 67 atoms the size error is of the order of 10 meV/atom. 

We tested for fc-point errors in 5Ut\i by performing calculations using both four and 32 Monkhorst-Pack Q| sampling 
points. Since explicit AIMD calculations with so many fc-points would be extremely expensive, we use the following 
procedure. From an existing F-point simulation we take a set of typically 10 atomic configurations separated by 0.1 ps. 
The ab initio total energies of these configurations calculated with the different /c-point samplings are then compared. 
For sampling with four fc-points, we did calculations on systems of up to 241 atoms, but the heavier calculations 
with 32 fc-points were done only on the 67-atom system. The results of these tests for the thermodynamic state 

V = 8.67 A-^/atom and T = 4300 K are also reported in table |l| where we see that for the smallest system containing 
67 atoms the difference with respect to a calculation with the F-point only is ~ 9 meV/atom, but as the number of 
atoms in the cell is increased above ^ 125 the difference becomes negligible. The result for the calculation with 32 
fc-points is identical to the one with four fc-points and is not reported in the table. We also found that fluctuations of 
the energy differences between the calculations done with the F-point only and those with 4 fc-points are extremely 
small. 

Similar, but less extensive, tests of system-size and fc-point errors have also been performed at the state V — 
6.97 A'^/atom, T = 6000 K, and we find that the variation of these errors with system size is numerically almost 
the same as before. The indication is therefore that 5Ut\i can be obtained to a precision of ca. 5 meV/atom from 
simulations on systems of 125 atoms or more. Unfortunately, it is not practicable yet to do all our AIMD simulations 
with this system size, and in practice we have computed 5Ut\i from F-point simulations on the 67-atom system, and 
corrected the results by adding 10 meV/atom, which from the present evidence appears to the almost constant error 
in the F-point 67-atoin results. 

As expected, the numerical values of 6Ut\i are small, and depend weakly on temperature and pressure across the 
range of thermodynamic states of interest. We find that they can be represented to within ^ 3 meV/atom by a sum 
of third-degree polynomials in V and T: 

3 

,5f/th/Af = ^ (a,F' + 6,r) , (19) 

with the following fitting parameters (units of eV, A and K): ag = 0.649; ai = —4.33 x 10^^; 02 = —4.19 x 10^^^; 
as = 6.48 X 10-5; 60 = 0.296; 61 = -6.51 x 10"^ 62 = 7.46 x 10"^; 63 = -2.07 x lO^^^. 



To test the validity of neglecting the higher-order fluctuation terms omitted from Eqn (15), we have performed 
full thermodynamic integration for four different thermodynamic states, the first three with V = 8.67 A^/atom and 
T = 4300, 6000 and 8000 K and the fourth with V = 6.97 AVatom and T = 8000 K, using the five equally spaced 
A values 0, 0.25, 0.5, 0.75 and 1.0. These calculations were done using F-point sampling on the system of 67 atoms. 
We have seen that this system size is not big enough to yield the required precision for Fai, but it should certainly 



be enough to test the adequacy of the second-order formula. In Table III we report a comparison between the results 
obtained from the integral using the five A values and those from the second order formula, and we see that they 
are practically indistinguishable. The Table also indicates that the term {{5/S.UY) ai /"^k-oT is rather insensitive to 
thermodynamic state and can be approximated to the required precision by setting it equal to 10 meV/atom. We 
have used this constant value in evaluating the ab initio free energy by Eqn dlq). 

IV. MELTING PROPERTIES 

From our parameterized formulas for the ab initio Helmholtz free energies F{V,T) of the h.c.p. solid (Ref. [ [ll[ ) 
and the liquid (present work), we immediately obtain the Gibbs free energies G{p,T) = F{V,T) +pV, and for each 



pressure the melting temperature T^ is determined as the T at which the latter free energies are equal for the solid 
and liquid. The resulting melting curve is reported in Fig. |l| for pressures from 50 to 350 GPa. On the same plot, 
we show the ab initio melting curve reported very recently by Laio et al. ||2^. We also compare with experimental 
melting curves or points obtained by shock experiments or by static-compression using the diamond anvil cell (DAC). 
DAC determinations of the melting curve of Fe and other transition metals have been performed by several research 
groups [|l7|-p0[. The early results of Williams et al. [^ lie considerably above those of other groups, and are now 
generally discounted. This still leaves a range of ca. 400 K in the experimental T,n at 100 GPa. Even allowing for 
this uncertainty, we acknowledge that our melting curve lies appreciably above the surviving DAC curves, with our 
Tm being above that of Shen et al. [|l9| by ca. 400 K at 100 GPa. We return to this discrepancy below. 

Shock measurements should in principle be able to fix a point on the high-pressure melting curve at the thermody- 
namic state where melting first occurs on the Hugoniot. However, temperature is notoriously difficult to measure in 
shock experiments. The temperatures obtained by Yoo et al. [ p2[ using pyrometric techniques are generally regarded 
as being too high by at least 1000 K. This has been confirmed by our recent ab initio calculations |ll| of Hugoniot 
temperature for h.c.p. Fe. We therefore disregard their data point on the melting curve. In the shock measurements 
of Brown and McQueen |2^ and Nguyen and Holmes [gj| , no attempt was made to measure temperature, which was 
estimated using models for the specific heat and Griineisen parameter; the approximate validity of these models is 
supported by our ab initio calculations ||ll[] on h.c.p. Fe. However, the identification of the Hugoniot melting point 
has been hampered by the possible existence of a solid-solid transition. In their measurements of sound velocity on 
the Hugoniot, Brown and McQueen ||2^ believed that they had observed a solid-solid transition as well as a separate 
melting transition. The new shock results of Nguyen and Holmes [M using improved techniques indicate that there 
is no solid-solid transition, and we place greater weight on their Hugoniot melting point. We plot in Fig. |l| the point 
reported by Brown and McQueen p3] as lying on the melting curve, though for the reasons just explained, we are 
cautious about accepting it. We also plot the point obtained from the measurements of Nguyen and Holmes |gj]. The 
pressure of 221 GPa is taken directly from their measurement of the onset of melting, while the temperature at this 
point is taken from our calculation of the Hugoniot temperature of the h.c.p. solid at this pressure, as reported in 
Ref. [^ (see also following section). 

We now consider possible sources of error in our DFT calculations. First, we recall that even with the best available 
GGA for exchange-correlation energy the low-temperature p{V) relation for h.c.p. Fe is not in perfect agreement with 



experiment. This has been shown by a number of independent calculations using all-electron techniques |30 3lf| as well 
as pseudopotential g^ and PAW [^,^ techniques, all of which agree closely with each other. Roughly speaking, the 
pressure is underpredicted by ca. 10 GPa at near-ambient pressures and by ca. 8 GPa in the region of 300 GPa. The 
pressure error can be thought of as arising from an error in the Helmholtz free energy, so that the true free energy 
Ftrue can be written as Ftruc — Fgga + SF, where -Fgga is our calculated free energy and 6F is the correction. If we 
take the pressure error 5p = —{dSF/dV)T to be linear in the volume, then 6F can be represented as 5F = biV + b2V^, 
where &i and 62 are adjustable parameters determined by least-squares fitting to the experimental pressure. If we 
now neglect the temperature dependence of SF, and simply add 6F{V) to the calculated free energies of solid and 
liquid, this gives a way of gauging our likely errors. We find that this free-energy correction leads to a lowering of the 
melting curve by ca. 350 K in the region of 50 GPa and by ca. 70 K in the region of 300 GPa. 

The second error source we consider is the PAW implementation, and specifically our choice of the division into 
core and valence states, and the PAW core radii. As mentioned earlier, at Earth's core pressures the 3p electrons, 
and to a lesser extent the 3s electrons, must be treated as valence states. Moreover, the choice of the maximum 
PAW core radius may also affect the calculations, because under such high pressures and temperatures the atoms 
come so close that the cores overlap. These errors may affect the melting curve if they fail to cancel between the 
liquid and the solid. To check both these possible problems, we have performed trial PAW calculations with the much 
smaller core radius of 0.85 A and with both 3s and 3p states in the valence set; with this choice of core radius the 
overlap of the cores in the liquid and the high temperature solid is almost negligible. We have then used Eq. (||) 
to calculate the free energy difference between the systems described with the two PAW approximations, repeating 
the calculations for both the liquid and the solid. To do that we have drawn two sets of 30 statistically independent 
configurations from two long simulations performed with the original PAW approximation on the solid and the liquid 
at V — 7.18 A'^/atom and T = 6700 K. As expected, we find a significant shift in the total electronic (free) energies. 
This shift is almost constant, thus validating the use of Eq. (H), but the important result is that it is almost the same 
for the liquid and the solid, the two numbers being F^^^^ — F^^^^. — —0.210 eV/atom F^aixi ~ ^soft = —0.204 eV/atom. 
Here, F^^^^ is the free energy calculated with small core and 3s and 3p states in valence, and F^^^^ the free energy 
with large core and the 3s and 3p frozen in the core, plus the effective pair-potential; the superscripts s and / indicate 
the solid and the liquid respectively. The effect is small, and stabilises the liquid by 6 meV/atom, which has the effect 
of shifting the melting curve down by '^ 60 K. 

As we show in Fig. H^, if we include both these corrections they bring our low-temperature melting curve into 
quite respectable agreement with the DAC measurements of Shen et al., while leaving the agreement with the shock 



points of Ng uyen and Holmes essentially unaffected. There is still a considerable discrepancy with the DAC curve of 



Boehler 17 1 and the ab initio results of Laic et al. |25| ] 

We now turn to the changes of volume and entropy on melting. Our calculated volume of melting (volume of liquid 
minus volume of coexisting h.c.p. solid at each pressure expressed as a percentage of the volume of the solid at that 
point) is plotted as a function of pressure in Fig. 0. We also show the melting volume predicted by the ab initio 
calculations of Laio et al. [g5| at the pressure 330 GPa, and it is encouraging to note that their value of 1.6 % is quite 
close to ours. The free-energy correction discussed above makes only a small difference to the calculated volume of 
melting: at 50 GPa the correction makes the volume of melting increase from 5.0 to 5.8 %, while at 300 GPa it is 
affected by less than 0.1 %. The most striking feature of our results is the steep decrease of Ay by a factor of about 
three in the range from 50 to 200 GPa, and its approximate constancy after that. 

Our predicted entropy of melting ASm (entropy per atom of liquid minus entropy per atom of coexisting solid) is 
plotted as a function of pressure in Fig. ||, where we also show the ab initio value of Laio et al. [^ at 330 GPa. The 
agreement of our value (1.05 fce) with theirs (0.86 fee) is reasonably close. The entropy of melting also decreases with 
increasing p, but more moderately than AV/V, the decrease between 50 and 200 GPa being only 30 %. We note 
the relevance to the slope of the melting curve, given by the Clausius-Clapeyron relation: dTm/dp = AV/AS. (This 
relation is satisfied identically by our results, since they are all derived from free energies.) The strong decrease of 
dTm/dp between 50 and 200 GPa and its approximate constancy thereafter is mainly due to the variation of AV/V. 

V. HUGONIOT PROPERTIES 

Since shock experiments are the only direct way of obtaining thermodynamic information for high-p/high-T liquid 
Fe, it is important to test our predictions against the available shock data. The data that emerge most directly from 
shock experiments consist of a relation between the pressure pn and the molar volume Vh on the so-called Hugoniot 
line, which is the set of thermodynamic states given by the Rankine-Hugoniot formula pSj : 

Iph{Vo - Vk) =En-Eo, (20) 

where En is the molar internal energy behind the shock front, and Eq and Vq are the molar internal energy and 
volume in the zero-pressure state ahead of the front. The pressure- volume and temperature-pressure relations on the 
Hugoniot are straightforwardly obtained from our ab initio calculations: for a given Vh, one seeks the temperature 
Th at which the Rankine-Hugoniot relation is satisfied, and from this one obtains pn (and, if required, E^i). In 
experiments on Fe, Vq and Eq refer to the zero-pressure b.c.c. crystal. We obtain Eq directly from GGA calculations 
that we performed on the ferromagnetic b.c.c. crystal, as described earlier [ |ll[ , but we use the experimental value 
of Vq. The slight shift produced by using instead the theoretical value of Vq was noted earlier |I^. Melting in shock 
experiments is usually detected by monitoring the sound velocity ||2^,^ , which shows marked discontinuities of slope 
along the Hugoniot. In a simple melting transition, there are discontinuities at two characteristic pressures ps and 
pi, which are the points where the solid and liquid Hugoniots meet the melting curve. Below ps, the material behind 
the shock front is entirely solid, while above pi it is entirely liquid; between ps and p;, the material is a two-phase 
mixture. 

We present in Fig. ^ our calculated Th(ph) Hugoniot curve for the liquid, together with our curve for the solid 



reported earlier [[111 and our ab initio melting curve. Without the free energy correction SF of Sec. IV, we find 
Ps = 229 and pi — 285 GPa. The very recent shock data of Nguyen and Holmes ||2J] give values of 221 and 260 GPa 
respectively, so that our ps value is very close to theirs, and our pi value is also not very different. 

If the correction SF is included in calculating the melting curve, then for consistency it must be included also in 
the solid and liquid Hugoniots. It is straightforward to obtain the corrected pn and £'h as a function of Vn for the 
two phases. But in the Rankine-Hugoniot equation we also need Eq for the b.c.c. crystal, and this will be subject to 
a correction similar to SF, but of unknown size. To supply the missing information, we add to the ab initio energy 
of b.c.c. Fe a correction term SF^cc, which we represent as ci -I- C2V. The constants ci and C2 are fixed by requiring 
that the equilibrium volume of the b.c.c. crystal and the low-temperature transition pressure between the b.c.c. and 
h.c.p. phases be correctly given. The resulting 'corrected' Th(ph) Hugoniots of the solid and liquid are reported in 
Fig. ^. The shifts in the curves are of about the same size as those discussed earlier [Q for the solid when we replaced 
the calculated b.c.c. volume Vq in the Rankine-Hugoniot equation by its experimental value, and are an indication 
of the inherent uncertainty due to DFT errors. The corrected values Ps = 243, pi = 298 GPa are now in somewhat 
poorer agreement with the experimental values. This is a rather sensitive test of DFT errors, since the shallow angle 
at which the Hugoniot curves cross the melting line amplifies the effect of the errors. 

We now turn to our hquid-state results for ph(Vh) compared with the shock data of Brown and McQueen |2^] 
(Fig. H), including for completeness our results for the solid reported earfier Q. We report results both with and 



without the free energy correction 6F, using the experimental b.c.c. volume Vq in the Rankine-Hugoniot equation in 
both cases. We mark on the Figure the volumes above which the shocked material is entirely solid and below which 
it is entirely liquid. Above the upper volume, we report our calculated h.c.p. Hugoniots, and below the lower volume 
the liquid Hugoniots. In the interval between them, we linearly mix the two. We note that the SF correction makes 
little difference to the liquid Hugoniot, which lies above the experimental values by ca. 3 %. 

Shock experiments on Fe have given values for the adiabatic sound speed vs — [Kg/ pY^"^ of the liquid, with Kg the 
adiabatic bulk modulus and p the mass density. Fig. o shows our ah initio values for vs of the liquid as a function of 
pressure on the Hugoniot, both with and without the 5F correction, compared with the shock data of Refs. |2^. Up 
to the pressure of ^ 260 GPa, the experimental points refer to the solid or the two-phase region, so it is in the liquid 
region above this pressure that the comparison is significant. In that region, our agreement with the experimental 
data is close, the discrepancies being ~ 2 and < 1 % for our uncorrected and corrected vs values respectively. 

We conclude this Section by reporting results for the Griineisen parameter 7 on the liquid Hugoniot. This parameter 
is defined as 7 = V{dp/dE)v — aKTV-mlC^^ with a the volume expansion coefficient, Kt the isothermal bulk 
modulus, C^ the constant-volume molar specific heat, and V^-a the molar volume. Assumptions or estimates of its 
values have played a key role in constructing parameterised equations of state for Fe. Our calculated 7 on the liquid 
Hugoniot is almost exactly constant, varying in the narrow range from 1.51 to 1.52 as p goes from 280 to 340 GPa. 

VI. THERMODYNAMICS OF THE LIQUID 

Although directly measured data on high-p/high-T liquid Fe all come from shock experiments, attempts have been 
made to combine these data with measurements at lower p and T using parameterised models for quantities such as 
Ks^ 7 and C« to estimate thermodynamic properties away from the Hugoniot curve [Q. These attempts have been 
crucial in trying to understand how the properties of the Earth's liquid core deviate from those of pure liquid Fe. We 
present here a brief comparison with these experimentally based extrapolations for the two quantities that determine 
the seismic properties of the outer core: the density p and the adiabatic bulk modulus Ks- 

Since the outer core is in a state of turbulent convection, the variation of its thermodynamic properties with depth 
is expected to follow an adiabat. We therefore present our comparisons on adiabats specified by their temperature 
TicB at p = 330 GPa, which is the pressure at the inncr-core/outer-core boundary (ICB) B^. We choose the two 
temperatures Ticb = 5000 and 7000 K, because the results of Sec. |^ indicate that the melting temperature at the 
ICB pressure lies between these limits. Our comparisons (Tables |^ and M) show that the uncorrected ah initio 
density is very close (within a few tenths of a percent) to the extrapolated experimental data at p = 150 GPa, and 
in slightly poorer agreement (within ~ 1.5 %) at p = 350 GPa. As expected, the free-energy correction lowers the 
predicted density, resulting in larger discrepancies with experiment of 1.5 % and 2.5 % at p = 150 and 350 GPa 
respectively. Our uncorrected ah initio Kg values also agree more closely with the experimental data, being typically 
within 2 %, while the corrected predictions disagree with the data by up to 8 %. However, given the closer agreement 
between ah initio and experiment on the Hugoniot (Sec. M), it is possible that some of the disagreements may be due 
to deficiencies in the experimental extrapolation. 

VII. DISCUSSION AND CONCLUSIONS 

In assessing the reliability of our results, we consider three sources of error: first, the uncontrolled DFT errors 
inherent in the GGA for exchange and correlation energy; second, the controllable errors in the detailed electronic- 
structure implementation of GGA, and specifically in the use of PAW to calculate the total ah initio (free) energy 
t/Ai(R.i, • ■ ■ R-Tv; T'oi) for each set of atomic positions Ri,...RAr; third, the statistical-mechanical errors, including 
system-size effects. We have endeavoured to reduce errors of the third kind below 10 meV/atom for the liquid. In 
our earlier free-energy calculations on the h.c.p. solid ||ll|], the corresponding error was estimated as ^ 15 meV/atom. 
Taking these errors together, and recalling that the resulting error in melting temperature Tm is roughly the combined 
free-energy error divided by Boltzmann's constant, we find an expected T^ error of ca. ±300 K. We have also attempted 
to control errors of the second kind by changing the division between core and valence states and by reducing the core 
radius. These tests suggest that the associated error in T^ is probably no more than ca. ±100 K. The inherent DFT 
errors are more difficult to quantify, but we have demonstrated that the known discrepancies in the low-temperature 
p(y) relation for h.c.p. Fe almost certainly lead to an overestimate of T,„ by ca. 350 K at 50 GPa and ca. 70 K at 
300 GPa, and we have corrected for this. We have also seen the significant shifts in the Hugoniot curves resulting 
from DFT errors. We believe the remaining uncertainty in T„j from this source could be as much as 300 K. 



Our attempts to correct for DFT errors give a melting curve which is in quite good agreement with the recent 
measurements of Shen et al. |19| , and with estimates based on shock data ||23|] ; the methods used to estimate tem- 
perature in the shock experiments are also supported by our ab initio results for Griineisen parameter 7 and specific 
heat Cv |l^] . Our melting curve is still above the experimental data of Boehler jl^ by ^ 800 K in the pressure region 
up to ca. 100 GPa. We cannot rule out the possibility that some of this discrepancy is due to our DFT errors. Our 
substantial disagreement with the ab initio melting curve of Laio et al. psj must be due to other reasons. We are 
currently working with authors of Rcf. [ p5[ to discover the cause of the disagreement, and we hope to report on this 
in the future. 

A key part of our strategy for eliminating system-size errors in the calculated free energies is the use of an empirical 
reference model which accurately reproduces the fluctuations of total energy. At first sight, the use of a reference 
model based on a purely repulsive pair potential might seem surprising, since it does not explicitly include a description 
of metallic bonding. An empirical reference model (there called an 'optimised potential model') is also used in the 
work of Laio et al. pq ], though they use it in a different way from us. Their optimised potential model is a form of 
the 'embedded atom model' (EAM) p8|-|5(|, which explicitly includes metallic bonding. As described in our earlier 
work |1^], we have investigated the consequences of using the EAM as a reference model. We showed there that for 
present purposes fluctuations of the bonding energy are negligible, and that under these circumstances the EAM is 
almost exactly equivalent to a model based on repulsive pair potentials. We also showed that there is no numerical 
advantage in using the EAM as reference model for the calculation of free energies. The use of different reference 
models per se therefore appears to have nothing to do with the current disagreement between ab initio melting curves. 

The agreement of our ah initio results with the limited data from shock experiments on the liquid is reasonably 
satisfactory. In particular, our predicted Hugoniot relation ph(^h) is almost as good as we found earlier for the solid. 
The adiabatic sound velocity of the liquid is also predicted to within 1 — 3 %, the discrepancy depending on whether 
or not we attempt to correct for DFT errors. The good agreement for the Griineisen parameter 7 is also encouraging. 
Our results for the h.c.p. solid ||ll| indicated that 7 varies little with pressure or temperature for 100 < p < 300 GPa 
and 4000 < T < 6000 K, and has a value of ca. 1.5. Our present results indicate that the same is true of the liquid. 

The ab initio free-energy techniques outlined here could clearly be adapted to a wide range of other problems, so 
that melting curves could be calculated for many materials, including those of geological interest, like silicates. We 
have recently completed ab initio calculations of the melting curve of aluminium up to pressures of 150 GPa, which 
are in excellent agreement with static-compression and shock data, as will be reported elsewhere pit]. 

In conclusion, we have shown how ab initio free-energy calculations based on thermodynamic integration can be 
used to obtain the melting curve and the volume and entropy of melting of a material over a wide pressure range. 
We have emphasised that the key requirement on the reference system used in thermodynamic integration is that it 
faithfully mimics the fluctuations of ab initio energy in thermal equilibrium. Our ab initio melting curve of Fe over 
the pressure range 50 — 350 GPa agrees fairly well with experimental data obtained from both static-compression and 
shock techniques, but significant discrepancies remain to be resolved. Our ab initio predictions for quantities obtained 
directly from shock experiments, including the Griineisen parameter of the liquid, agree closely with the measured 
data in most cases. 
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P (kg m-3) 


r(K) 


9540 


10700 


11010 
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13300 
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0.097 










4300 


(60) 


0.085 








5000 




(132) 
0.089 








6000 


0.104 


(140) 
0.096 


0.089 


0.103 


0.125 


7000 


(90) 


(151) 
0.093 


(170) 
0.098 


(251) 
0.109 


(360) 
0.131 


8000 




(161) 
0.092 


(181) 
0.099 


(264) 
0.104 


(375) 
0.124 






(172) 


(191) 


(275) 


(390) 



TABLE I. Normalised fluctuation strength a (see text) characterising the accuracy with which the inverse-power reference 
model mimics the energy fluctuations of ah initio hquid Fe. Values of a (eV units) are given for a set of AIMD simulations at 
different densities p and temperatures T. Pressure at each thermodynamic state (GPa units) is given in parenthesis. 





(<5(7th(iV)-5C/th(241))/iV(cV) 


(<5(7tt 


(4fe) - 5f/th(r))/7V (eV) 


a (eV) 


67 


-0.009 ± 0.002 




0.009 ±0.003 


0.085 


89 


-0.012 ± 0.001 




0.007 ±0.002 


0.073 


107 


-0.010 ± 0.001 




0.006 ±0.002 


0.083 


127 


0.004 ± 0.001 




0.000 ±0.002 


0.086 


157 


0.001 ± 0.001 




0.001 ±0.002 


0.069 


199 


0.001 ± 0.001 






0.082 


241 


0.000 ± 0.001 




0.001 ±0.002 


0.101 



TABLE II. Dependence on number of atoms A'' in the simulation ceU of size errors and fc-point samphng errors in the quantity 
5[/th entering the ah initio free energy of hquid Fe (see Eqn (llq)). Second column reports 5Ut\i/N (eV units) with a constant 
offset chosen so that the reported value for the largest system size is zero. Third column reports difference of SUth/N between 
simulations using four fc-points and F-point sampling. Fourth column reports normalised fluctuation strength a (see text) for 
different system sizes. 



r(K) 



N-^ ll d\ (AC/);, (eV) 



{{5AUf)Ai/2NkBT (eV) 



4300 
6000 
8000 



0.012 

0.010 

0.006 (0.010) 



0.012 

0.009 

0.006 (0.010) 



TABLE III. Difference _Fai — ^ref between free energies of ah initio and reference systems calculated in two ways: by full 
thermodynamic integration (column 2) as in Eqn (kh, and by the second-order fluctuation approximation (column 3) as in 
Eqn (o). Free energy differences are given per atom in eV units for three temperatures at the density p = 10700 kg m~^. Values 
in parenthesis refer to the density p = 13300 kg m"''. 
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12180 


300 


12844 


(12756) 




12575 


(12481) 




13000 






12800 


350 


13315 


(13232) 




13043 


(12970) 




13550 






13290 



TABLE IV. Comparison of ab initio and experimental density p of liquid Fe on two adiabats, with adiabats specified by the 
temperature T at the pressure p = 300 GPa. Ab initio p values are given both without and with (in parenthesis) free-energy 
correction SF (see text). 
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P (GPa) 


T = 5000 K 




T = 7000 K 




T = 5000 K 
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= 7000 K 


150 


708 (662) 




656 (613) 




695 




668 


200 


878 (838) 




820 (781) 




877 




849 


250 


1050 (1010) 




981 (944) 




1058 




1016 


300 


1220 (1180) 




1140 (1103) 




1232 




1193 


350 


1384 (1350) 




1296 (1264) 




1400 




1355 



TABLE V. Comparison of ab initio and experimental adiabatic bulk modulus Ks of hquid Fo on two adiabats, with adiabats 
specified by the temperature T at the pressure p = 330 GPa. Ab initio Ks values are given both without and with (in 
parenthesis) free-energy correction 5F (see text). 
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FIG. 1. Comparison of melting curve of Fe from present calculations with previous experimental and ab initio results: heavy 
solid and dashed curves: present work without and with free-energy correction (see text); chain curve: ab initio results of 
Ref. psl]; dots, light dashes and squares: DAC measurements of Refs. pif , ||l7| and [M; triangles, diamond and solid square: 
shock experiments of Refs. p3], p3| and |24]. Error bars are those quoted in original references. 
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FIG. 2. Ab initio fractional volume change on melting of Fe as a function of pressure. Solid and dashed curves: present 
work, without and with free-energy correction (see text); black dot: Ref. |E3]. 
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FIG. 3. Ab initio entropy change on melting per atom (units of Boltzmann's constant fee). Solid and dashed curves: present 
work, without and with free-energy correction (see text); black dot: Ref. feq]. 
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FIG. 4. Relation between ah initio melting curve and ah initio Hugoniot temperature-pressure curves. Heavy continuous 
and dashed curves: melting curves calculated without and with free-energy correction (see text); light continuous and chain 
curves: Hugoniot of solid without and with free-energy correction; light dashed and dotted curves: Hugoniot of liquid without 
and with free-energy correction. 
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FIG. 5. Ab initio Hugoniot pressure- volume curve compared with experimental results of Ref. p3|. Solid and dashed curves: 
ab initio results without and with free-energy correction (see text); squares: experimental results. Vertical dotted lines indicate 
volumes at which melting starts and finishes according to present (uncorrected) ab initio results. To the right of rightmost 
vertical dotted line, curves represent solid Hugoniot from Ref. IllJ; to the left of leftmost vertical line, curves represent present 
liquid Hugoniot. 
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FIG. 6. Longitudinal speed of sound on the Hugoniot. Circles: experimental values from Ref. 
curves: present ab initio values without and with free-energy correction (see text). 
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